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Abstract 

The  U.S.  Army  Engineer  Research  and  Development  Center  (ERDC), 
Coastal  and  Hydraulics  Laboratory  (CHL),  has  undertaken  the 
development  of  the  multimodule  Adaptive  Hydraulics  (ADH) 
hydrodynamic,  sediment,  water  quality,  and  transport  numerical  code. 
The  Mellor-Yamada  level  2.5  and  k-s  turbulence  closure  models  were 
incorporated  into  ADH-SW3.  This  report  documents  the  incorporation  of 
these  turbulence  closure  models  into  ADH-SW3;  it  also  documents  the 
validation  of  their  incorporation  by  using  them  in  an  application  of  the 
ADH-SW3  model  code  to  three  flume  experiments.  These  flume  tests  were 
designed  to  ensure  that  the  incorporated  turbulence  closure  models  in 
ADH-SW3  are  solving  the  pertinent  equations  accurately.  A  basic 
description  of  other  turbulence  closure  options  implemented  with  ADH- 
SW3  is  also  provided,  and  a  user  manual  (User  Appendix)  is  included. 
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Preface 

This  report  presents  the  details  of  incorporating  the  Mellor-Yamada  level 
2.5  and  k-8  turbulence  closure  schemes  into  the  three-dimensional  shallow 
water  module  of  the  Adaptive  Hydraulics  (ADH-SW3)  numerical  code.  The 
report  also  briefly  describes  other  turbulence  models  incorporated  into 
ADH-SW3  and  includes  a  user  manual  (User  Appendix).  This  investigation 
was  conducted  from  October  2013  through  March  2014  at  the  U.S.  Army 
Engineer  Research  and  Development  Center  (ERDC)  by  Dr.  Gaurav  Savant 
of  Dynamic  Solutions  EEC. 

This  report  is  published  as  a  product  of  the  Elood  and  Coastal  Storm 
Damage  Reduction  Program  of  the  ERDC,  Vicksburg,  MS.  Dr.  Cary  Talbot 
was  the  Program  Manager,  and  William  Curtis  was  the  Technical  Director. 

At  the  time  of  publication  of  this  report.  Dr.  Jeffery  P.  Holland  was 
Director  of  ERDC,  and  ETC  John  T.  Tucker  HI  was  Acting  Commander. 
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Unit  Conversion  Factors 


Multiply 

By 

To  Obtain 

cubic  feet 

0.02831685 

cubic  meters 

cubic  yards 

0.7645549 
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0.01745329 
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feet 
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meters 
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0.5144444 

meters  per  second 

microns 

1.0  E-06 
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meters 

miies  (U.S.  statute) 

1,609.347 

meters 

miies  per  hour 

0.44704 

meters  per  second 

pounds  (force) 

4.448222 

newtons 

siugs 

14.59390 

kiiograms 

square  feet 

0.09290304 

square  meters 

square  miies 

2.589998  E+06 

square  meters 

square  yards 

0.8361274 

square  meters 

yards 

0.9144 
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1  Introduction 

The  U.S.  Army  Corps  of  Engineers  (USAGE),  through  the  U.S.  Army 
Engineer  Research  and  Development  Center  (ERDC),  has  developed  a 
robust  multidimensional  mass  conservative  finite  element  hydrodynamic 
and  constituent  transport  numerical  code  Adaptive  Hydraulics  (ADH). 
Adaptive  Hydraulics  has  been  referred  to  as  “ADH”  and  “AdH”  in 
literature;  the  abbreviation  ADH  is  used  in  this  report  in  accordance  with 
how  Adaptive  Hydraulics  is  referenced  in  peer-reviewed  literature. 

The  objectives  of  this  study  were  to  incorporate  the  Mellor-Yamada  level 
2.5  (MY-25)  and  k-8  turbulence  closure  schemes  into  ADH-SW3,  document 
the  incorporation  and  the  verification  and  validation  of  the  schemes,  and 
provide  a  user  manual  (User  Appendix).  Both  MY-2.5  and  k-8  models 
involve  the  solution  of  two  transport  equations  to  represent  the  generation 
and  dissipation  of  turbulent  energy.  These  two  models  have  the  capability  to 
represent  several  real-world  phenomena  such  as  density  stratification  and 
adverse  pressure  flow.  Additionally,  other  closures  exist  in  the  ADH-SW3 
code  —  the  Mellor-Yamada  level  2  and  the  Smagorinsky  —  and  a  brief 
description  of  those  is  presented  for  comparison  and  completeness. 
Validation  of  these  closures  was  documented  in  Savant  et  al.  (2014). 
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2  Equation  Development 

Basic  equations 

A  variety  of  three-dimensional  (3D)  reservoir,  riverine,  coastal,  and 
estuarine  simulation  models  (e.g.,  ADH-SW3,  CH3D,  TABS-MDS)  exist. 
These  models  solve  some  form  of  the  Reynolds-averaged  Navier-Stokes 
equations  using  the  hydrostatic  assumption.  The  governing  equations  for 
velocity  17,  =  [U,V,W) ,  and  constituent  C  in  x,.  —  (x,y,z)  is  written  as 


dU. 

- ^  =  0 

dXj 


dU-  d(up.)  1  dP 

1 

1 

f  d 

dt  dXj  Pq  ^x^ 

Po 

^dz 

m  o  i 


■F.  =0 


wherej  is  summed  over  (x,  y,  z)  or  (U,  V,  W)  and  i  is  (x,  y)  or  (U,  V),  as 
appropriate;  the  equation  for  the  z-direction  component  reduces  to 


■I 

P{z)  =  Pa+  f 


(1) 

(2) 

(3) 


under  the  hydrostatic  pressure  assumption,  where  is  the  atmospheric 
pressure,  g  is  the  acceleration  due  to  gravity,  p(z)  is  the  density  at  a 
location  in  the  z-direction,  and  q  is  the  water  surface  elevation. 


In  addition,  the  convection-diffusion  equation  for  a  baroclinic  transport 
constituent  is  written  as 


dC 

dt 


dz  ^  dz 


,^+F^=0 


(4) 


with  the  sum  over  j  as  before. 

Coefficients  for  the  horizontal  viscosity  [F^ )  and  horizontal  diffusivity  (F^) 
terms  are  generally  calculated  using  a  horizontal  closure  scheme. 
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Smagorinsky  for  ADH-SW3.  The  vertical  eddy  viscosity  for  momentum 
{K^=E^  =  )  and  eddy  diffusivity  for  transport  iK^=D^)  coefficients 

are  utilized  to  represent  the  Reynolds  stresses  and  turbulence  transport. 
Estimates  of  these  coefficients  are  obtained  by  turbulence  closure  schemes, 
which  are  discussed  in  the  following  sections. 

Horizontal  eddy  viscosity  and  diffusivity  -  Smagorinsky 

Coefficients  for  the  horizontal  eddy  viscosity  (Ej  )  and  horizontal  diffusivity 
( Ep )  terms  are  generally  calculated  using  a  horizontal  closure  scheme. 


These  terms  are  usually  represented  as 


Pi=}- 

Po 


_d_E„^ 
dX;  dX; 


dx.  ^dx. 


(5) 

(6) 


wherej  is  summed  over  (x,  y)  or  (C7,  V)  and  i  is  (x,  y),  as  appropriate.  ADH- 
SW3  applies  the  Smagorinsky  (1963)  technique  to  calculate  the  coefficients 
Pxx  j  Pxy  ’  Pyy  ’  Pyx  >  Dy .  Thc  formulation  takes  the  form 


E  =E  =E  =D  =D  =cA 

XT  yy  xy  x  y  '■'s  , 


[9t7' 

2 

_L 

(dV^ 

2 

_L 

'dU' 

2 

_L 

[dV'' 

1 

2 

+  2 

'dUdV  dUdV] 

dx , 

^dx , 

\ 

.dy , 

dy  dx  dx  dy 

(7) 


where  is  the  Smagorinsky  coefficient  (o.05<  <o.i)  and  A  is  the 

element  area. 


Generic  length  scale  (GLS)  implementation 

To  facilitate  the  incorporation  of  additional  turbulence  closure  schemes, 
MY-2.5  and  k-8  in  ADH-SW3  were  implemented  using  the  GLS  scheme 
proposed  by  Umlauf  and  Burchard  (2003).  The  GLS  scheme  is  a  two- 
equation  turbulence  model  that  takes  advantage  of  the  similarities 
between  most  two-equation  models.  The  first  equation  represents  the 
transport  and  evolution  of  the  standard  turbulent  energy  (k),  and  the 
second  equation  represents  the  transport  and  evolution  of  a  generic 
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quantity  (\|/).  A  generic  quantity  (\|/)  is  used  to  establish  the  turbulence 
length  scale  (/).  The  equation  for  k  is  represented  as 


dk 

dt 


+  U, 


dk 

d 

\K^ 

dk 

dXj 

dz 

dz 

/ 

-(P  +  B-^)  =  0 


(8) 


where j  is  summed  over  (x,  y,  z)  or  (U,  V,  W),  as  appropriate,  o,^is  the 

turbulence  Schmidt  number  for  k  (Table  i),  P  and  B  represent  the 
generation  of  turbulence  by  shear  and  buoyancy  production  or  destruction 
by  density  gradients,  respectively,  and  8  is  the  dissipation  of  turbulence.  P, 
B  and  8  are  represented  as 


with 


P  =  K^M^ 


M"  = 


dU 


dz 


+ 


dV 


\dz  , 


B  =  -K,N^  , 

Po  9z 


^  p  3  m  -1 

3+^  ,  -H - - 


8  =  (c°)  "ip" 


(9) 

(10) 

(11) 

(12) 


where  Po  is  the  density  at  standard  temperature  and  pressure  (STP),  c°  is  a 

model  stability  coefficient,  and  m,  n,  andp  are  model  specihc  parameters 
(Table  i). 

The  second  equation  describes  the  evolution  and  transport  of  a  generic 
parameter  (\|/),  and  is  specihed  as 


+u,  - 

d 

\K^dw] 

dt 

^dx. 

dz 

[a,  9z] 

-^(c,P  +  C35-C28F^„„)  =  0  (13) 


wherej  is  summed  over  (x,  y,  z)  or  (U,  V,  W),  as  appropriate,  and  C3 
are  model  specific  parameters  (Table  1),  o^  is  the  turbulence  Schmidt 
number  for  \\f  (Table  1),  and 
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=  (14) 

The  turbulent  length  scale  ( / )  is  recovered  through 

/  =  (c°f^V  (15) 

The  MY-2.5  model  requires  the  specification  of  a  wall  function  ), 

several  formulations  of  which  are  presented  in  the  literature.  Details  of  the 
wall  function  are  provided  later  in  this  report. 

Careful  selection  of  parameters  p,  m,  n,  Fwaii,  Ci,  C2,  C3,  o^.,  and  results 

in  the  recovery  of  exact  formulations  for  the  MY-2.5  and  k-s  turbulence 
closure  schemes. 


Table  1.  GLS  Parameters 


k-kl 

k-e 

p 

0.0 

3.0 

m 

1.0 

1.5 

n 

1.0 

-1.0 

2.44 

1.0 

2.44 

1.3 

Cl 

0.9 

1.44 

C2 

0.5 

1.92 

C3 

0.9 

1.0 

0 

0.5544 

0.5544 

Mellor-Yamada  level  2.5  (IVIY-2.5)  GLS  implementation 

The  MY-25  scheme  is  reproduced  from  (8)  and  (13)  using  the  k-kl  column 
in  Table  1  and  specifying  p  =  0.0,  m  =  1.0  and  n  =  1.0,  yielding 


dt  ' 


dk 

d 

dx- 

dz 

dk 

dz 


-(P  +  B-e)  =  0 


(16) 


and 
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dkl  ^  dkl  d 
dt  ^  dx-  dz 


K, 


dkl 


dz 


-  +  C3E  -  C28F^,„  )  =  0 


The  original  implementation  of  MY-2. 5  defined  as 


Ko=>S„M 


where  =  0.2. 

In  the  GLS  implementation  of  MY-25,  is  defined  as 


(17) 


(18) 


i^n=—  (19) 

where  is  the  Schmidt  number  for  k  and  takes  a  value  of  2.44,  Ci 

=  0.9,  C2  =  0.5,  and  C3  =  0.9. 

The  wall  function  can  be  prescribed  as  originally  suggested  in  Mellor- 
Yamada  (1982) 

Fu^,  =  2h-d,  (20) 


or  as 


77  — 

-^wall  — 


1  + 


KMIN(d^,d,) 


(21) 


as  suggested  by  Burchard  et  al.  (1998)  or  as 


77  — 


1  +  Ej 


11. 

Kd 


(22) 


as  suggested  by  Burchard  (2001)  for  deep  basins,  or  as  the  following  which 
was  suggested  by  Blumberg  et  al.  (1992)  for  general  open  channel  flows 
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77  — 

^wall  ~ 


1  +  E2 

\  ^  ] 

2 

+  E, 

/  f] 

(23) 


In  the  ADH-SW3  implementation,  each  of  the  wall  functions  is  imple¬ 
mented,  where  =  1.33,  =  0.25,  k  is  the  von  Karman  constant  and 

takes  a  value  of  0.4,  h  is  the  flow  depth,  is  distance  to  surface,  and  is 

distance  to  the  bed.  Experience  has  shown  that  Equation  (23)  is  the  wall 
function  that  is  likely  most  useful  for  estuaries,  and  Equation  (21)  is  more 
appropriate  to  riverine  applications. 

The  vertical  eddy  viscosity  for  momentum  ( )  and  eddy  diffusivity  for 
transport  ( )  coefficients  are  obtained  by  appealing  to  formulations 
involving  stability  functions.  These  are  represented  as 


K„  =  SJ^ 

(24) 

K,  =  S,lM 

(25) 

where  Sm  and  Sh  are  stability  functions  for  momentum  and  transport, 
respectively,  and  are  derived  algebraically  (Warner  et  al.  2005). 

There  are  several  formulations  of  Sm  and  Sh  available,  such  as  Kantha  and 
Clayson  (1994),  Galperin  et  al.  (1988),  and  Canuto  et  al.  (2001)  among 
several  others,  in  the  literature.  The  ADH-SW3  implementation  of  MY-25 
uses  the  formulation  described  by  Kantha  and  Clayson  (1994)  and  is 
described  as 


Sn  = 


6A, 


B. 


1  J 


I-3A2G;,  6\+BJl-C. 


(26) 


5.=  A 


l-^-3C, 

5, 


9\2A,+AJl-C,)]S,G, 


1-9A,A,G, 


(27) 


with 
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(/)  2 

G,^=^N^  ,  -  0.28  <G;,  <0.0233  (28) 

k 

Kantha  and  Clayson  (1994)  stability  function  parameter  values  are 
presented  in  Table  2. 


Table  2.  MY-25  stabilization  function  parameters. 


Parameter 

4 

4 

4 

4 

4 

Q 

G,, 

min 

G,, 

max 

Value 

0.92 

0.74 

16.6 

10.1 

0.08 

0.7 

0.2 

-0.28 

0.0233 

The  GLS  formulation  imposes  an  upper  limit  on  the  length  scale  (/)  to 
account  for  the  reduction  in  mixing  in  stable  stratification.  Warner  et  al. 
(2005)  present  this  limit  as 

p^O^  iV">0  (29) 

-  at" 

Vertical  eddy  viscosity  through  Mellor-Yamada  level  2  (IVIY-2) 

ADH-SW3  has  the  option  of  computing  the  coefficients  and  using 

the  MY-2  scheme  described  in  Mellor  and  Yamada  (1982).  This  closure 
scheme  is  a  simplification  of  the  MY-2.5  scheme:  it  neglects  diffusion, 
equates  the  generation  of  turbulent  energy  with  its  dissipation,  and  assumes 
horizontal  homogeneity.  In  essence,  it  is  a  turbulence  equilibrium  model 
(i.e.,  any  turbulence  generated  or  transported  is  instantaneously 
dissipated).  This  implies  that 


P  +  B  =  £  (30) 

Equation  (30)  and  Mellor  and  Yamada’s  (1982)  assumption  that  turbulent 
energy  is  destroyed  as  it  is  created,  and  thus  not  transported,  in  essence 
removes  Equations  16  and  17  from  the  computation.  (Eor  a  detailed 
discussion  of  the  MY-2  scheme,  see  Mellor  and  Yamada  (1982)). 

MY-2  does  not  follow  the  assumptions  in  the  GLS  scheme,  and  to  obtain 
the  MY-2  formulation,  the  basic  equations  developed  by  Mellor  and 
Yamada  (1982)  must  be  used.  These  equations  include  defining  the  length 
scale  (/)  as 
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l  =  Kd, 


(31) 


or  as 


/ 


Kd, 


(32) 


where  all  variables  are  as  previously  defined.  Equation  31  (Mellor  and 
Blumberg  2004)  is  widely  used  and  provides  a  linear  shape  for  the  mixing 
length.  ADH-SW3  uses  equation  32  similar  to  that  suggested  by  Robert 
and  Ouellet  (1987).  In  the  absence  of  surface  waves,  Equation  32  goes  to  o 
at  the  surface  and  at  the  bed.  If  the  master  length  scale  is  assumed  to  be 
specified  as  above  or  by  some  other  relationship,  a  flux  Richardson 
number  (R/)  is  defined  as 


ml 


(33) 


The  gradient  Richardson  number  (R)  is  defined  as 


R  = 


~M^ 


or 


R 


R 


f 


(34) 


(35) 


Mellor-Yamada  (1982)  defined  the  following  as  the  energy  budget  for  the 
turbulent  kinetic  energy: 


This  can  be  further  written  as 


(36) 


1-R. 


(37) 
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Utilizing  the  above  and  definitions  for  and  ,  the  turbulent 
momentum  and  constituent  transport  coefficients  are  written  as 


and 


=  SP 


du 

dz 

(!-«/)] 

0.5 


(38) 


du 

r  /  \i 

Kn  =  Sf 

dz 

0.5 


(39) 


K-e  GLS  implementation 

The  k-8  model  is  retrieved  from  the  GLS  model  described  by  Equations  (8) 
through  (15)  by  utilizing  the  parameters  specified  in  the  k-s  column  of 
Table  1.  The  equations  for  turbulent  kinetic  energy  (k)  and  the  dissipation 
(£■)  are  written  in  GLS  as 


dt  ^ 


dk  d 

\K^dk] 

1 

[o^  dz] 

-(P  +  B-z)  =  0 


(40) 


where  all  quantities  are  the  same  as  previously  described.  ol  is  the 
Schmidt  number  for  eddy  diffusivity  of  k  and  has  a  value  of  1.0,  and 


—+u. 

dt  ^ 


dz 

d 

de 

dXj 

dz 

dz 

/ 

■— (c.P  +  CoE-c^e)  =  0 

k 


(41) 


where  is  calculated  using  Equations  12, 14, 15  and  24,  and  o^,  =  o^^,  is 
the  eddy  diffusivity  for  8  with  a  value  of  1.3. 

Boundary  specification 

Two-equation  turbulence  models  do  not  resolve  the  viscous  sublayer,  so 
boundary  conditions  are  applied  in  order  for  the  models  to  accurately 
represent  the  turbulence  transport  processes.  The  boundary  conditions  for 
the  MY-25  and  k-8  two-equation  turbulence  models  can  be  specified  by 
assuming  that  production  of  turbulence  by  shear  (P)  is  completely  balanced 
by  the  dissipation  of  that  turbulence  (£■).  Generation  of  turbulent  kinetic 
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energy  through  buoyancy  is  negligible  since  density  gradients  tend  towards 
zero. 

The  resulting  boundary  conditions  for  k  are  (Warner  et  al.  2005) 


where  k  is  the  turbulent  kinetic  energy  as  before,  c°  is  the  model  specific 
constant  (Table  1),  u  is  the  shear  velocity,  and  the  subscripts  b  and  s 
represent  bed  and  surface,  respectively. 

The  resulting  boundary  conditions  for  T  are  (Warner  et  al.  2005) 


where  p,  m,  and  n  are  model  specific  constants  (Table  1),  k  is  the  von 
Karman  constant,  and  and  are  the  mixing  depths  for  bottom  and 

surface,  respectively.  Stability  concerns  require  that  the  value  of  is 
specified  as  0.01  (in  the  units  of  model),  and  z^  is  specified  as  1%  of  the 
depth.  The  shear  velocities  for  the  surface  and  bed  are  calculated  using 


7^1^  Z„ 

In - 1  +  -^ 


,  and  u*  = 


/n  —  -1  +  ^ 


where  u  is  the  depth  averaged  flow  velocity,  is  the  wind  velocity,  h  is 
the  flow  depth  ,  Zq  is  the  bottom  roughness  height,  and  z^g  is  the  surface 
roughness  height  and  is  equal  to  Zg  in  the  absence  of  surface  waves  and  to 

0.85  times  the  root  mean  square  (RMS)  wave  height  in  the  presence  of 
surface  waves  (Mellor  and  Blumberg  2004). 

Minimum  vaiues  of  turbuience  generation  and  dissipation 

Turbulence  generation  and  dissipation  values  appear  in  the  denominator 
terms  for  several  computed  quantities;  this  necessitates  the  specification 
of  minimum  values  for  both.  Experience  has  shown  that  these  minimum 
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values  are  application  specific,  and  the  user  must  determine  the  minimum 
values  required  for  stability. 

Buoyancy  suppression  functions 

It  is  well  known  that  stable  stratification  reduces  momentum  transfer 
across  a  pycnocline,  and  several  formulations  exist  to  suppress  the  eddy 
viscosity  and  diffusivity  across  this  high-density  gradient  region.  ADH- 
SW3  provides  a  choice  of  five  functions  to  suppress  eddy  viscosity  and 
diffusivity,  each  defined  as  a  modification  of  and  . 


Henderson-Sellers 

Henderson  and  Sellers  (1982)  derived  a  suppression  formulation  based  on 
observations  of  the  atmospheric  boundary  layer,  and  is  represented  as 


K„ 


1  +  0.1  AR 


(45) 


1  +  31 


(46) 


Munk-Anderson 

Munk  and  Anderson  (1948)  were  some  of  the  first  researchers  to  observe 
and  detail  the  influence  of  stable  stratification  on  momentum  transfer 
across  a  pycnocline  based  on  observations  of  oceanic  thermoclines.  They 
provided  the  following  modifier  to  the  vertical  eddy  viscosity  and  diffusivity: 


^  _ _ yvn _ 

[1  +  P(n0i?r 

K 

'  [1  +  P(n,)i?p 


(47) 


(48) 


Values  of  parameters  Ui,  112,  (3(n^)  and  (3(n2)  are  provided  in  Table  3. 
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Table  3.  Buoyancy  suppression  function  parameters. 


ni 

na 

ADH-SW3 

Option 

Henderson-Sellers 

1 

1 

0.74 

37 

1 

Munk-Anderson 

0.5 

1.5 

10 

3.3 

2 

Kent-Pritchard 

2 

- 

0.24 

- 

3 

Pritchard 

2 

- 

0.28 

- 

4 

French- 

McCutcheon 

2 

" 

10 

- 

5 

Kent-Pritchard 

Kent  and  Pritchard  (1957)  derived  their  suppression  function  from  data 
collected  in  the  James  River  Estuary  and  is  represented  as 


K, 


K. 


m,h 


m,h 


l  +  P(n^)R 


(49) 


Pritchard 

Pritchard  (i960)  re-evaluated  observations  from  the  James  River  Estuary 
and  represented  his  formulation  as 


K, 


K. 


m,h 


m,h 


(50) 


French-McCutcheon 

Erench  and  McCutcheon  (1983)  derived  their  suppression  formulation 
from  observed  data  in  the  Great  Ourse  Estuary,  and  their  formulation 
takes  the  form 


K, 


K. 


m,h 


m,h 


l  +  P(ni)R 


(51) 


Note  that  though  Equations  (47)  to  (51)  have  similar  forms,  the  values  of 
|3(n)  vary  by  more  than  an  order  of  magnitude.  Eurthermore,  these 

functions  are  based  on  empirical  data,  and  care  should  be  exercised  in 
their  application. 
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3  Testing 

Testing  of  MY-2.5  and  k-s  was  performed  on  three  test  cases.  The  first  case 
tested  the  generation/dissipation  of  turbulence  due  to  high-velocity 
gradients  in  the  absence  of  density  gradients;  the  second  case  tested  the 
generation/dissipation  of  turbulence  primarily  due  to  high  density 
gradients;  the  third  case  tested  the  performance  of  the  model  for  flow  past 
a  cylindrical  obstruction  in  the  flow  path.  The  third  test  case  is  a 
particularly  difficult  test  for  any  hydrostatic  code  because  of  the  vertical 
accelerations  involved;  in  the  absence  of  robust  turbulence  models,  the 
code  will  fail  to  simulate,  in  particular,  the  bottom  currents  or  the 
downstream  vortices.  MY-2  and  Smagorinsky  closure  schemes  have  been 
tested  previously  (Savant  et  al.  2014)  and  will  not  be  considered. 

Flow  around  an  emergent  spur  dike 

This  test  case,  based  upon  the  work  presented  in  Rajaratnam  and 
Nwachukwu  (1983),  is  designed  to  test  the  accuracy  and  adequacy  of  the 
turbulence  closure  schemes  implemented  into  ADH-SW3. 

The  test  domain  is  illustrated  in  Figure  1.  An  emergent  spur  of  0.152 
meters  (m)  length  and  0.03  m  width  is  placed  14.0  m  downstream  of  the 
inflow  location  (at  the  left  boundary).  A  uniform  flow  of  0.0453  cubic 
meters  per  second  (m3/s)  is  applied  at  the  left  boundary  with  a  tail  water 
elevation  of  0.189  m  applied  at  the  right  boundary. 

Figure  1.  Domain  for  spur  dike  test. 


The  model  parameters  utilized  are  as  follows: 

Smagorinsky  coefficient:  0.2 

Uniform  background  eddy  viscosity:  0.0015  square  meters  per  second  (m  /s) 
Manning’s  n  value:  0.01. 

Minimum  k:  0.000005  m^s'^ 
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Minimum  \|/  or  s:  0.0005  m^s'^ 

Number  of  vertical  layers;  4. 

Figures  2  and  3  show  the  recirculation  zone  computed  by  MY-25  and  k-8, 
respectively.  The  recirculation  lengths  computed  were  11.8  and  12.3  by 
MY-25  and  k-s,  respectively.  These  values  match  the  value  of  12  times  the 
spur  length  reported  in  literature  (Wang  et  al.  2009). 

Figure  2.  Recirculation  zone  and  velocity  magnitudes  with  MY-25  (the  color  ramps  from  slow 
to  fast  velocity,  with  red  representing  velocities  >  0.25  m/s). 


Figure  3.  Recirculation  zone  and  velocity  magnitudes  with  k-s  (the  color  ramps  from  slow  to 
fast  velocity,  with  red  representing  velocities  >  0.25  m/s). 


The  average  velocity  in  the  domain  at  steady  state  was  computed  by  the 
model  to  be  0.257  m/s  using  MY-25  and  0.282  m/s  using  k-s.  A  close 
match  between  the  two  provides  further  confidence  that  both  MY-25  and 
k-s  were  implemented  correctly  for  a  high-velocity  gradient  case. 

Propagation  of  salinity  subsequent  to  a  lock  exchange 

This  test  was  run  to  ascertain  the  ability  of  the  model  to  accurately  represent 
the  speed  (CT)  of  a  density  wedge,  referred  to  as  the  shock  speed  in  Shin  et  al. 
(2004).  The  test  consisted  of  a  2  m  long,  0.2  m  wide,  and  0.2  m  deep  flume 
with  denser  salt  water,  35  parts  per  thousand  (ppt),  in  the  left  half  and 
freshwater,  o  ppt,  in  the  right  half.  The  barrier  separating  the  two  is 
instantaneously  removed  allowing  the  denser  fluid  to  slump  under  the 
lighter  fluid  and  move  as  a  density  wedge.  As  in  Shin  et  al.  (2004),  U  is 
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determined  by  noting  the  time  (f)  for  the  salinity  to  increase  a  certain 
amount  a  distance  (x),  from  the  initial  separating  barrier:  U  =  x/t.  Figures  4 
and  5  illustrate  the  domain  and  initial  constituent  state,  respectively,  for 
this  test. 


Figure  4.  Domain  for  iock  exchange. 


Figure  5.  Initiai  constituent  state  for  lock  exchange. 


The  model-computed  shock  speed  is  used  to  calculate  the  densiometric 
Froude  number  as 


79(1^ 


(52) 


where  y  is  the  ratio  of  lower  density  to  higher  density  (0.997  for  this  test) 
and  h  is  the  total  dense  fluid  depth.  The  iy,  computed  for  this  test  case  is 

0.51  for  MY-25  and  0.5  for  k-s;  Shin  et  al.  (2004)  reported  the  value  of  0.5 
for  the  energy-conserving  value  of  nonrigid  lid  density  currents  (Shin  et  al 
2004),  as  calculated  with  ADH-SW3  here. 
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The  model  parameters  utilized  are  as  follows: 

Smagorinsky  coefficient:  0.2 

Uniform  background  eddy  viscosity:  0.0000001  mVs 
Manning’s  n  value:  0.01 
Minimum  k:  0.000005  ni^s  ^ 

Minimum  \|/  or  s:  0.00000001  m^s  ^ 

Buoyancy  suppression  function:  Henderson-Sellers 
Number  of  vertical  layers:  5 

Figures  6  and  7  illustrate  the  results,  showing  the  position  of  the  red 
salinity  wedge,  for  the  MY-25  and  k-s  closure  schemes  tests  after  time  10  s. 
It  is  observed  that  MY-25  is  slightly  more  diffusive  compared  to  k-s.  The 
reason  for  this  diffusivity  is  unclear  at  present  and  is  a  subject  of 
investigation. 


Figure  6.  State  of  density  wedge  at  10  s  with  MY-25  scheme. 


Figure  7.  State  of  density  wedge  at  10  s  with  k-s  scheme. 


Flow  around  a  cylindrical  obstruction 

This  test  case,  roughly  based  upon  the  work  of  Melville  and  Raudkivi 
(1977),  was  designed  to  test  the  adequacy  of  the  turbulence  closure 
schemes  to  qualitatively  model  flow  conditions  where  flow  separation 
occurs  and  vertical  accelerations  are  significant.  ADH-SW3  is  a  hydrostatic 
code  and  assumes  negligible  vertical  accelerations;  hence,  the  burden  of 
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accounting  for  these  vertical  accelerations  rests  entirely  on  the  turbulence 
closure  schemes  implemented. 

Figure  8  illustrates  the  domain  for  this  test  case.  A  solid  cylinder,  5  m  in 
diameter,  is  placed  at  a  distance  of  200  m  from  the  upstream  end,  and 
flow  has  to  go  around  the  cylinder  causing  flow  separation  and  plunging  of 
flow.  The  inflow  is  specified  as  2000  m3/s,  and  a  constant  tail  water 
elevation  of  10  m  is  specified  at  the  downstream  end. 


Figure  8.  Domain  for  flow  around  a  cyiindertest. 


The  model  parameters  utilized  are  as  follows: 

Smagorinsky  coefficient:  0.2 

Uniform  background  eddy  viscosity:  0.0000001  mVs 

Manning’s  n  value:  0.01 

Minimum  k  for  MY-25:  o.  oooooim^s  ^ 

Minimum  \\f  for  MY-25  or  s:  0.5  m^s  s 

Minimum  k  for  k-s:  o.  05  m^s  ^ 

Minimum  \|/  or  s  for  k-s:  0.05  m^s  s 

Buoyancy  suppression  function:  Kent-Pritchard 
Number  of  Vertical  Layers:  5 

Melville  and  Raudkivi  (1977)  present  results  from  mobile  bed  experiments 
and  present  some  velocity  results  from  a  case  before  initiation  of  bed  scour. 
The  observed  results  show  that  flow  reaches  stagnation  upstream  of  the 
cylinder,  and  downstream  flow  separation  and  reattachment  occurs  (on  the 
cylinder,  opposite  of  the  upstream  stagnation  point).  Figures  9  and  10 
provide  the  bottom  streamtraces  for  simulated  results  from  the  MY-25  and 
k-8  schemes,  respectively.  Comparing  the  simulation  results  with  those 
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provided  by  Melville  and  Raudkivi  (1977,  Figure  2,  where  velocity  results 
from  a  case  before  initiation  of  bed  scour  [i.e.,  a  flat  bed]  are  shown),  it  is 
noted  that  both  schemes  adequately  represent  the  bottom  velocity  behavior. 
The  k-s  scheme  provides  a  better  representation  of  the  downstream  eddies, 
whereas  the  MY-25  scheme  is  too  diffusive  and  results  in  smaller/ weaker 
eddies.  In  the  absence  of  exact  experimental  data,  it  is  not  possible  to 
comment  quantitatively  on  the  performance  of  the  two  schemes. 

Figure  9.  Bottom  streamtrace,  MY-25.  Red  indicates  higher  veiocity  and  biue  iower. 


Flow 


ERDC/CHL  CR-15-1 


20 


Figure  10.  Bottom  streamtrace,  k-e.  Red  indicates  higher  veiocity  and  biue  lower. 
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4  Summary  and  Conclusions 

This  report  presents  the  development  and  incorporation  of  generalized 
length  scale  based  Mellor-Yamada  level  2.5  (MY-25)  and  k-8  turbulence 
closure  schemes  into  the  3D  shallow  water  module  of  ADH  (ADH-SW3). 
Also  presented  are  the  details  of  other  turbulence  schemes  available  within 
ADH-SW3  (Mellor  Yamada  level  2  and  Smagorinsky),  the  wall  functions 
implemented  for  the  MY-25  scheme,  the  buoyancy  suppression  functions 
available,  and  details  of  three  test  cases. 

The  three  test  cases  were  designed  to  exercise  the  turbulence  closure 
schemes  under  a  high-velocity  gradient,  under  a  high-density  gradient, 
and  under  significant  vertical  accelerations.  Both  MY-25  and  k-8  were  able 
to  accurately  represent  the  physics  of  the  problems.  It  was  observed  that 
for  the  high-density  gradient  test  case,  the  MY-25  was  slightly  more 
diffusive  than  the  k-8  closure  scheme.  The  reason  for  this  is  under  active 
investigation;  the  preliminary  thinking  is  that  the  difference  in  diffusivity 
is  due  to  the  wall  function  utilized  in  MY-25(the  (k-e)  scheme  does  not 
require  a  wall  function). 

Future  work  in  turbulence  closure  is  required,  and  it  is  suggested  that  an 
option  be  included  in  ADH-SW3  to  write  out  nodal  eddy  viscosity  and 
diffusivity  values.  The  capability  to  write  out  the  model-computed  eddy 
viscosity  values  by  node  is  essential  to  visualize  model  behavior;  this  is 
especially  important  to  perform  an  analysis  of  density  stratification.  This  is 
not  possible  at  present  because  ADH-SW3  considers  eddy  viscosity  and 
diffusivity  as  elemental  properties. 
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Appendix;  User’s  Guide  for  Turbuience 
Options 

Turbulence  cards 

The  turbulence  options  in  ADH-SW3  are  activated  using  the  “MP  TUR” 
card  and  the  turbulence  dirichlet  boundaries  are  specified  using  “DB  TRI” 
cards.  Specification  of  boundary  conditions  for  turbulence  constituents  is 
not  mandatory  but  is  recommended  for  stability  purposes  and  consistency 
with  published  literature. 

If  a  turbulence  scheme  other  than  Mellor-Yamada  level  2  (MY-2)  is  used, 
the  user  must  specify  “CN  TKE”  and  “CN  TDS”  for  the  turbulent  kinetic 
energy  and  turbulence  dissipation,  respectively.  The  format  of  these  cards 
is  as  follows: 

CN  TKE  “Transport  number”  “Reference  concentration” 

CN  TKE  “Transport  number”  “Reference  concentration” 

Turbulence  control  cards 

MPTUR 

TURBULENCE  PARAMETERS 


Field 

Type 

Value 

Description 

1 

char 

MP 

Card  type 

2 

char 

TUR 

Parameter 

3 

int 

>  0 

Material  number  (not  string) 

4 

int 

>  0 

Horizontal  turbulence 

5 

int 

any  int 

Vertical  turbulence 

6 

real 

any  real 

Horizontal  turbulence  coefficient 

7 

int 

>0 

Buoyancy  suppression  function 
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If  the  vertical  turbulence  option  specified  is  Mellor-Yamada  2.5  or  option 
“2”,  an  additional  card  for  the  wall  function  must  be  specified  and  is 

8  int  1-4  Wall  proximity  function 

If  the  vertical  turbulence  option  specified  is  greater  than  “1,”  additional 
card  values  must  be  specified  and  are 


9 

real 

>0 

Minimum  turbulent  kinetic  energy  (TKE) 

10 

real 

>0 

Minimum  turbulent  dissipation  (TDS) 

10 

real 

>0 

Maximum  mixing  length  fraction 

Table  Ai  presents  field  4,  5,  7,  and  8  value  options  that  activate  various 
turbulence  schemes  as  well  as  the  buoyancy  suppression  functions  in 
ADH-SW3. 


Table  A 1.  Field  parameter  values. 


Field  4 

Field  5 

Field  7 

Field  8 

Smagorinsky  -  0 

Mellor-Yamada  Level 
2-1 

Henderson-Sellers  - 1 

Mellor-Yamada  -  1 

Water  Surface 

Slope  - 1 

Mellor-Yamada  Level 
2.5-2 

Munk-Anderson  -  2 

Burchard  et  al.  -  2 

k-e-3 

Kent-Pritchard  -  3 

Burchard  -  3 

No  vertical 
turbulence  <  0 

Pritchard  -  4 

Blumberg  et  al.  -  4 

French-McCutcheon  -  5 

Since  turbulence  options  are  material  region  specific  parameters,  the  user 
can  specify  different  turbulence  options  and/or  parameters  for  various 
material  regions  within  the  simulation. 
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DB  TRI 

Field 

Type 

Value 

TURBULENCE  BOUNDARY 

Description 

1 

char 

DB 

Card  type 

2 

char 

TRI 

Parameter 

3 

int 

>  0 

Node  string  number 

4 

int 

>  0 

Transport  constituent 

5 

int 

0  or  1 

Surface  or  wall/bottom 

6 

int 

any int 

Vertical  turbulence  option 

Guidance  for  buoyancy  suppression  function 

The  buoyancy  suppression  function  utilized  has  a  direct  and  substantial 
impact  on  the  diffusion  of  momentum  and  energy  across  the  pycnocline; 
this  is  especially  true  for  stable  stratification.  The  choice  of  which  function 
to  use  is  based  on  the  type  of  environment  being  modeled. 

Testing  performed  using  ADH-SW3  has  provided  the  following  basic 
guidance  for  selection  of  the  buoyancy  suppression  function: 

1.  Henderson-Sellers  is  best  suited  for  systems  that  do  not  have  stable 
stratification  or  for  a  relatively  mixed  estuarine/riverine  system.  Examples 
include  estuaries  such  as  Galveston  Bay. 

2.  Munk-Anderson  and  French-McCutcheon  are  suited  for  systems  that  have 
a  substantial  freshwater  input  into  a  large  estuary  and  have  density 
stratification.  Examples  include  estuaries  such  as  Mobile  Bay. 

3.  Kent-Pritchard  and  Pritchard  are  not  recommended  except  for 
experimental  testing  due  to  a  lack  of  confidence  in  the  data  utilized  to 
develop  the  buoyancy  suppression  function. 

Since  turbulence  options  are  material  region  specific  parameters,  the  user 
can  specify  different  turbulence  options  and/or  parameters  for  various 
material  regions  within  the  simulation. 
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